clc;clear all;close all;
global alpha betta depr sigmaC sigmaL phi govexp omega taumax N T ALPHA popweight
global kstst lstst cstst rstst wstst kstart kistart tauK0
global Delta1 Delta2 lamda transfer kg lam kmax
global restrict

%run mainbench.m;

load('benchN2.mat')
%load('benchN5.mat')

kistart = kj;
kstart = sum(popweight.*kistart)+kg;
kmax=3;

lamda=lam;
restrict = 0
Delta1 = .3;
Delta2 = Delta1;

kstst = fzero('findstst',kstart+0.3,options)

rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end

T=150; %we assume that the economy reaches the steady state after T periods

tauK0 = taumax;

%initial guesses 
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
    lb(i)=kstart-0.2;
    ub(i)=kstst+0.3;
    X(T+i) = lstst; %l
    lb(T+i)=0.15;
    ub(T+i)=0.4;
    if (N==5)
    X(2*T+i) = max(2-.07*i,0); %gamma
    end
    if (N==2)
    X(2*T+i) = max(1-.05*i,0); %gamma
    end
    lb(2*T+i)=0;
    ub(2*T+i)=5;
end;
transfer=1;
X(3*T+1:3*T+1+2*(N-1)) = [Delta1 transfer lamda];
lb(3*T+1:3*T+1+2*(N-1)) = [-1 -1*ones(1,N-1) 0.2*ones(1,N-1)];
ub(3*T+1:3*T+1+2*(N-1)) = [1 1*ones(1,N-1) 1.2*ones(1,N-1)];

ALPHall=[1 0.95 0.9 0.8 0.7 0.6 0.5 0.488 0.487 0.4865 0.4863 0.4862 0.4861 0.48 0.46 0.43 0.4 0.37 0.35 0.347 0.3467 0.3466 0.34 0.33 0.327 0.325 0.32 0.31 0.3095 0.309 0.3085 0.308 0.306 0.304 0.302 0.3]
ii=1;    
%capitalist max U at 0.3467 (ii=21)
%worker max U at 0.4861 (ii=13)

for ii=1:length(ALPHall)
clc;clear all;close all;
load('solN2_fb.mat')
ii=36
ALPHA=ALPHall(ii)
    
    if (ii>1)
        X=XX_fb(ii-1,:);
        if ii==8 | ii==9 | ii==19 | ii==31
            X=XX_fb(ii-2,:);
        end
        if ii==11 | ii==18 | ii==20 | ii==25
            X=XX_fb(ii-3,:);
        end
        if ii==21 | ii==28 | ii==29
            X=XX_fb(ii-4,:);
        end
        if ii==22 | ii==23 | ii==24 | ii==30 | ii==31
            X=XX_fb(ii-5,:);
        end
        if ii==32 | ii==33 | ii==34 | ii==35
            X=XX_fb(ii-8,:);
        end
        if ii==12 | ii==13 | ii==26
            X=0.5*XX_fb(ii-1,:)+0.5*XX_fb(ii-4,:);
        end
        if ii==36
            X=0.5*XX_fb(ii-6,:)+0.5*XX_fb(ii-12,:);
        end
        if ii==27
            X=0.5*XX_fb(ii-3,:)+0.5*XX_fb(ii-6,:);
        end
        
        Delta1 = X(3*T+1);
        Delta2 = Delta1;
        transfer = X(3*T+2:3*T+2+N-2);
        lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
    end
    for i=1:T;
        time1 = i/T;
        X(i) = time1*(kstst) + (1-time1)*kstart; %k
    end
kstst = fzero('findstst',kstart+0.3)
rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end

tauK0 = taumax;
tolf0 = 1e-7;

if ii==20
    [X,check] = broydn('resid3v',X,tolf0)
    if size(X,1)>1
        X=X';
    end
end

options = optimset('MaxFunEvals',10000000,'MaxIter',200,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
[X,check] = fsolve('resid3v',X,options)

if ii==1 | (ii>=9 & ii<=17) | ii==26
    X1=X;
    X=X1;
    for i=1:T;
        time1 = i/T;
        X(i) = time1*(kstst) + (1-time1)*kstart; %k
    end

    [X,check] = broydn('resid3v',X,tolf0)
    if size(X,1)>1
        X=X';
    end
    sum(resid3v(X).^2)
end

if ii==7 | ii==8
    [X,check] = broydn('resid3v',X,tolf0)
    if size(X,1)>1
        X=X';
    end
end

if ii==15 | ii==16 | ii==26
    options = optimset('MaxFunEvals',10000000,'MaxIter',200,'TolFun',1e-15,'TolX',1e-15,'Disp','iter');
    [X,check] = fsolve('resid3v',X,options)
end

sumresid_all(ii)=sum(resid3v(X).^2)

XX_fb(ii,:)=X;

%for ii=1:length(ALPHall)
%X=XX(ii,:);
k = X(1:T);
l = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
Delta2 = Delta1
transfer = X(3*T+2:3*T+2+N-2);
lamda = X(3*T+3+N-2:3*T+3+2*(N-2));
Delta1*kistart(1)+sum(Delta2.*kistart(2:N))
kstst = fzero('findstst',[kstart-0.2 kmax],options)%,0,0)

rstst = betta^(-1)-1+depr; %r from Euler at the steady state
estst = (rstst/alpha)^(1/(1-alpha))*kstst; % express total effective labor from r=F_k
for jj=2:N
ljpart(jj-1)=popweight(jj)*phi(jj)*lam(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL);
end
lstst=estst/(popweight(1)*phi(1)+sum(ljpart)); %express l1 using the expression for total effective labor and l2=f(lam,l1)
for jj=2:N
    L2a(jj-1) = lamda(jj-1)^(-sigmaC/sigmaL)*(phi(jj)/phi(1))^(1/sigmaL); %last terms in (11) and f_l
end
wstst = (1-alpha)*kstst^alpha*estst^(-alpha); %w=F_e
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst.^alpha*estst.^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst = omega*lstst^sigmaL*(1+sum(ALPHA.*lamda.^(-sigmaC).*phi(2:N)'/phi(1).*L2a)+(Delta1+sum(Delta2.*phi(2:N)'/phi(1).*L2a))*(1+sigmaL))/(wstst*(popweight(1)*phi(1)+sum(ljpart)));

klast = [kstart X(1:T-1)]; %k_{t-1}
e=[l' (l'*(lamda.^(-sigmaC/sigmaL).*(phi(2:N)/phi(1))'.^(1/sigmaL)))]*(popweight'.*phi);
e=e';

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];
lnext = [X(T+2:2*T) lstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
listst = [lstst, L2(lamda,lstst,phi(2:N))']; %hours worked of groups in the long run
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt), L2(lamda,l(tt),phi(2:N))']; %hours worked of groups during the transition
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)- omega*li.^(1+sigmaL)/(1+sigmaL)))...
    + betta^(T)/(1-betta)*(log(cistst) - omega*listst.^(1+sigmaL)/(1+sigmaL));
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC) - omega*li.^(1+sigmaL)/(1+sigmaL)))...
    + betta^(T)/(1-betta)*((cistst.^(1-sigmaC)-1)/(1-sigmaC) - omega*listst.^(1+sigmaL)/(1+sigmaL));
end

if sum(V>=Vsq)==2
    PI(ii)=1
else
    PI(ii)=0
end

residall(ii)=sum(resid3v(X).^2)
Vall(ii,:)=V
Vsq

%tax rates
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
plot(tauKt)
PI(ii)
tauKt
lamda_all(ii)=lamda;
Delta1_all(ii)=Delta1;
Delta2_all(ii)=Delta2;
term_all(ii)=1+ALPHA*lamda^(1-sigmaC)+(Delta1+Delta2*lamda)*(1-sigmaC);
kappa=lamda^(-sigmaC/sigmaL)*(phi(2)/phi(1))^(1/sigmaL);
terml_all(ii)=1+ALPHA*kappa^(1+sigmaL)+(Delta1+phi(2)/phi(1)*kappa*Delta2)*(1+sigmaL);
cst_all(ii)=cstst;
gammast_all(ii)=gamma(T);
must_all(ii)=mustst;

%foc for labor
klam=lamda^(-sigmaC/sigmaL)*(phi(2)/phi(1))^(1/sigmaL);
lamderl(ii)=-lamda^(-sigmaC/sigmaL)*(phi(2)/phi(1))^((1+sigmaL)/sigmaL)*betta^T*(1+sigmaL)*lstst^(sigmaL)/(sum(betta.^[0:T-1].*c.^(-sigmaC).*c)+betta^T/(1-betta)*cstst^(-sigmaC)*cstst...
    -sigmaC/sigmaL*lamda^(-(sigmaC+sigmaL)/sigmaL)*(phi(2)/phi(1))^((1+sigmaL)/sigmaL)*(sum(betta.^[0:T-1].*(-omega*l.^(sigmaL).*l))+betta^T/(1-betta)*(-omega)*lstst^sigmaL*lstst));
%foc for consumption
lamderc(ii)=-lamda*betta^T*(1-sigmaC)*cstst^(-sigmaC)/(sum(betta.^[0:T-1].*c.^(-sigmaC).*c)+betta^T/(1-betta)*cstst^(-sigmaC)*cstst...
    -sigmaC/sigmaL*lamda^(-(sigmaC+sigmaL)/sigmaL)*(phi(2)/phi(1))^((1+sigmaL)/sigmaL)*(sum(betta.^[0:T-1].*(-omega*l.^(sigmaL).*l))+betta^T/(1-betta)*(-omega)*lstst^sigmaL*lstst));
%end

save('solN2_fb.mat')

end

Vcap_fb=Vall(:,1);
Vwor_fb=Vall(:,2);
save('solN2_fb.mat')

